Authors:Martin Ostrowski, Deepa Varkey Version: 20180515
An introduction to R for plotting data on maps.
This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.
plot(cars)
cars is a pre-loaded dataset for demonstration purposes. To see what the dataset looks like execute ‘cars’ in a chunk
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.
cars
speed dist
1 4 2
2 4 10
3 7 4
4 7 22
5 8 16
6 9 10
7 10 18
8 10 26
9 10 34
10 11 17
11 11 28
12 12 14
13 12 20
14 12 24
15 12 28
16 13 26
17 13 34
18 13 34
19 13 46
20 14 26
21 14 36
22 14 60
23 14 80
24 15 20
25 15 26
26 15 54
27 16 32
28 16 40
29 17 32
30 17 40
31 17 50
32 18 42
33 18 56
34 18 76
35 18 84
36 19 36
37 19 46
38 19 68
39 20 32
40 20 48
41 20 52
42 20 56
43 20 64
44 22 66
45 23 54
46 24 70
47 24 92
48 24 93
49 24 120
50 25 85
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).
The aims of this practical are to develop core data handling, plotting and mapping skills using a programming lanaguage. R is the language that we will be using today. Other languages, such as ArcGIS, Matlab and Python, are also very capable of achieving the same, or better, results. Oh, by the way, R is free. Ultimately the choice of language will be yours. These notes are written in R Markdown format which will allow you to execute the code inline and produce a web friendly (.html) format preview.
The objectives of todays practical are to:
setup
(20-30 min)
getting started
maps (45 min)
exercises: practice customising plots
extending
(30 min)
(10 min)
To plot data on a map you first need a map to work from. We can find the information by loading two packages, maps, and mapdata. Information in this section was sourced from Jeff Bowman
install.packages('maps')
install.packages('mapdata')
library(maps)
library(mapdata)
Let’s try plotting a map
map() # low resolution map of the world
let’s limit the map region to Australia by addin limits for longitude x and latitude y using xlim=c(110,157), ylim=c(-48, -10)
map('worldHires', xlim=c(110,157), ylim=c(-48, -10), fill=T, col="grey")
map('worldHires', xlim=c(110,180), ylim=c(-48, -10), fill=T, col="grey", asp=1)
Now we can add data points.
Australia operates an Integrated Marine Observing System which served data from ocean observations through the Australian Ocean Data Network
load data from the Integrated Marine Observing System National Reference Stations. IMOS operates 7 National Reference Stations around the coast of Australia.The coordinates for each station are given in the data file NRSlatlon.csv
Before loading the data we need to make sure that we are in the right directory with the setwd() command
setwd("~/Desktop/R4Oceanography/")
nrs.sites<-read.table("NRSlatlong.csv", h=T, sep=",")
check the data
nrs.sites
NRS code lat long
1 Port Hacking PHB -34.1192 151.2267
2 Maria Island MAI -42.5967 148.2333
3 North Stradbroke Island NSI -27.3450 153.5620
4 Rottnest Island ROT -32.0000 115.4167
5 Yongala YON -19.3085 147.6184
6 Kangaroo Island KAI -35.8322 136.4473
7 Darwin DAR -12.4000 130.7680
Now try adding the data to the map using points() to add a plotting character pch= to an x,y coordinate
map('worldHires', xlim=c(110,180), ylim=c(-48, -10), fill=T, col="grey", asp=1)
points(nrs.sites$long, nrs.sites$lat)
These points are hard to see. We can make them bigger using the cex= option, and use a different plotting character (pch=16)
map('worldHires', xlim=c(110,180), ylim=c(-48, -10), fill=T, col="grey", asp=1)
points(nrs.sites$long, nrs.sites$lat, pch=16, cex=3)
This looks great.
Q how can we tell which site is which? Use text() to add text to an x,y position on a map
map ('worldHires', xlim=c(110,180), ylim=c(-48, -10), fill=T, col="grey", asp=1)
points (nrs.sites$long, nrs.sites$lat, pch=16, cex=5)
text (nrs.sites$long, nrs.sites$lat, nrs.sites$code, col="white")
setwd("~/Desktop/R4Oceanography/")
The working directory was changed to /Users/mostrowski/Desktop/R4Oceanography inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.
in16.sites<-read.table("IN16latlong.csv", h=T, sep=",")
Check the data
in16.sites
CTD type lat long col
1 IN16v04.001 Test cast -34.33 151.76 red
2 IN16v04.002 Eddy Transect -34.17 152.57 red
3 IN16v04.003 Eddy Transect -34.27 152.33 red
4 IN16v04.004 Eddy Transect -34.32 152.06 red
5 IN16v04.005 Eddy Transect -34.40 151.80 red
6 IN16v04.008 Process Stn Coastal -33.39 152.03 red
7 IN16v04.009 Process Stn Coastal -33.39 152.02 red
8 IN16v04.010 -33.01 152.46 blue
9 IN16v04.011 Port Hacking -34.13 151.24 black
10 IN16v04.012 Upwelling? -36.32 150.28 red
11 IN16v04.013 Process Stn- Coastal -36.25 150.30 red
12 IN16v04.014 Process Stn- Coastal -36.26 150.29 red
13 IN16v04.015 Transect - Southern -36.23 151.27 red
14 IN16v04.016 Transect - Southern -36.25 151.93 red
15 IN16v04.017 Transect - Southern -36.25 153.26 red
16 IN16v04.018 Transect - Southern -36.25 154.25 red
17 IN16v04.019 Transect - Southern -37.00 154.00 red
18 IN16v04.020 Not the right water -37.00 152.99 red
19 IN16v04.021 Process Stn - Tasman -36.99 153.89 red
20 IN16v04.022 Process Stn - Tasman -36.99 153.89 red
21 IN16v04.023 Deep Cast - Tasman -36.89 152.55 purple
22 IN16v04.024 -36.49 150.29 red
23 IN16v04.025 -34.84 151.34 red
24 IN16v04.026 Process Stn - EAC -30.56 153.67 red
25 IN16v04.027 Process Stn - EAC -30.64 153.62 red
26 IN16v04.028 Drift Study - EAC -30.81 153.44 orange
27 IN16v04.029 Drift Study - EAC -30.97 153.39 orange
28 IN16v04.030 Drift Study - EAC -31.51 153.33 orange
29 IN16v04.031 Drift Study - EAC -31.65 153.32 orange
30 IN16v04.031 Drift Study - EAC -31.64 153.32 orange
31 IN16v04.032 Drift Study - EAC -31.84 153.32 orange
32 IN16v04.033 Drift Study - EAC -31.98 153.30 orange
33 IN16v04.034 Process Stn - EAC -32.46 153.18 red
34 IN16v04.035 Process Stn - EAC -32.48 153.16 red
35 IN16v04.035 Process Stn - EAC -32.47 153.17 red
36 IN16v04.036 Drift Study - EAC -32.57 153.12 orange
37 IN16v04.037 Drift Study - EAC -32.59 153.07 orange
38 IN16v04.038 Drift Study - EAC -32.76 152.99 orange
39 IN16v04.039 Drift Study - EAC -33.00 152.91 orange
40 IN16v04.040 Drift Study - EAC -33.05 152.88 orange
41 IN16v04.041 Process Stn - Tasman -32.47 154.49 red
42 IN16v04.042 Process Stn - Tasman -32.48 154.48 red
43 IN16v04.043 Transect - Middle -32.47 154.10 red
44 IN16v04.044 Transect - Middle -32.46 153.70 red
45 IN16v04.045 Transect - Middle -32.46 153.40 red
46 IN16v04.046 Transect - Middle -32.29 152.88 red
47 IN16v04.047 Process Stn - EAC -29.15 154.48 red
48 IN16v04.048 Process Stn - EAC -29.16 154.48 red
49 IN16v04.049 Transect - Northern -28.00 153.78 red
50 IN16v04.050 Transect - Northern -28.02 154.03 red
51 IN16v04.051 Transect - Northern -28.01 154.35 red
52 IN16v04.052 Transect - Northern -28.00 154.78 red
53 IN16v04.053 Transect - Northern -28.00 155.03 red
54 IN16v04.054 Deep Cast -27.92 155.15 red
55 IN16v04.055 Process Stn - EAC -28.00 155.04 red
Q: How many stations are there?
nrow(in16.sites)
[1] 55
Plot the stations on a map. Hint: change the longitude range.
map ('worldHires', xlim=c(145,160), ylim=c(-42, -20), fill=T, col="grey", asp=1)
title(xlab = 'Longitude',
ylab = 'Latitude')
box()
grid()
points(in16.sites$long, in16.sites$lat, pch=16, col="red")
col= (“red”, “blue”, “black”)pch= (1-25)map ('worldHires', xlim=c(145,160), ylim=c(-42, -20), fill=T, col="grey", asp=1)
points(in16.sites$long, in16.sites$lat, pch=16, col=in16.sites$col, cex=2)
Ocean Data Products specific to Australia are hosted on the AODN website. However, for this practical we will browse data from the NASA Ocean Color website, via the level 3 data browser. From the pulldown menus select Terra MODIS 4µ nightime, 4km resolution, 8 day composite. When hovering your mouse over the desired date you will see 3 download options in the bottom left, middle and right corners. Click on te leftee corner to download the .nc file and move it to your R4Oceanography directory on the desktop.
The .nc file extension signifies it is a netcdf file Network Common Data Form. Netcdf4 format is commonly used for high density arrays of data, such as those used in oceanography, climatology, taxonomy and GIS applications.
To enable R to work with netcdf files additional libraries need to be added
install.packages('ncdf4')
library(ncdf4)
install.packages('oce')
library(oce) #A Package for Oceanographic Analysis
What does a netcdf file look like?
sst01<-nc_open("~/Desktop/R4Oceanography/T20162572016264.L3m_8D_SST4_sst4_4km.nc")
sst01
File ~/Desktop/R4Oceanography/T20162572016264.L3m_8D_SST4_sst4_4km.nc (NC_FORMAT_NETCDF4):
3 variables (excluding dimension variables):
short sst4[lon,lat] (Chunking: [1729,40]) (Compression: level 4)
long_name: 4um Sea Surface Temperature
scale_factor: 0.00499999988824129
add_offset: 0
units: degree_C
standard_name: sea_surface_temperature
_FillValue: -32767
valid_min: -1000
valid_max: 10000
display_scale: linear
display_min: -2
display_max: 45
byte qual_sst4[lon,lat] (Chunking: [1729,40]) (Compression: level 4)
long_name: Quality Levels, Sea Surface Temperature
_FillValue: -1
valid_min: 0
valid_max: 5
unsigned byte palette[eightbitcolor,rgb] (Contiguous storage)
4 dimensions:
lat Size:4320
long_name: Latitude
units: degree_north
_FillValue: -999
valid_min: -90
valid_max: 90
lon Size:8640
long_name: Longitude
units: degree_east
_FillValue: -999
valid_min: -180
valid_max: 180
rgb Size:3
eightbitcolor Size:256
65 global attributes:
product_name: T20162572016264.L3m_8D_SST4_sst4_4km.nc
instrument: MODIS
title: HMODIST Level-3 Standard Mapped Image
project: Ocean Biology Processing Group (NASA/GSFC/OBPG)
platform: Terra
temporal_range: 8-day
processing_version: 2014.0
date_created: 2016-10-06T10:24:39.000Z
history: l3mapgen par=T20162572016264.L3m_8D_SST4_sst4_4km.nc.param
l2_flag_names: LAND,~HISOLZEN
time_coverage_start: 2016-09-12T09:00:10.000Z
time_coverage_end: 2016-09-20T11:45:08.000Z
start_orbit_number: 89021
end_orbit_number: 89139
map_projection: Equidistant Cylindrical
latitude_units: degrees_north
longitude_units: degrees_east
northernmost_latitude: 90
southernmost_latitude: -90
westernmost_longitude: -180
easternmost_longitude: 180
geospatial_lat_max: 90
geospatial_lat_min: -90
geospatial_lon_max: 180
geospatial_lon_min: -180
grid_mapping_name: latitude_longitude
latitude_step: 0.0416666679084301
longitude_step: 0.0416666679084301
sw_point_latitude: -89.9791641235352
sw_point_longitude: -179.97917175293
geospatial_lon_resolution: 4.63831233978271
geospatial_lat_resolution: 4.63831233978271
geospatial_lat_units: degrees_north
geospatial_lon_units: degrees_east
spatialResolution: 4.64 km
number_of_lines: 4320
number_of_columns: 8640
measure: Mean
suggested_image_scaling_minimum: -2
suggested_image_scaling_maximum: 45
suggested_image_scaling_type: LINEAR
suggested_image_scaling_applied: No
_lastModified: 2016-10-06T10:24:39.000Z
Conventions: CF-1.6
institution: NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group
standard_name_vocabulary: NetCDF Climate and Forecast (CF) Metadata Convention
Metadata_Conventions: Unidata Dataset Discovery v1.0
naming_authority: gov.nasa.gsfc.sci.oceandata
id: T20162572016264.L3b_8D_SST4.nc/L3/T20162572016264.L3b_8D_SST4.nc
license: http://science.nasa.gov/earth-science/earth-science-data/data-information-policy/
creator_name: NASA/GSFC/OBPG
publisher_name: NASA/GSFC/OBPG
creator_email: data@oceancolor.gsfc.nasa.gov
publisher_email: data@oceancolor.gsfc.nasa.gov
creator_url: http://oceandata.sci.gsfc.nasa.gov
publisher_url: http://oceandata.sci.gsfc.nasa.gov
processing_level: L3 Mapped
cdm_data_type: grid
identifier_product_doi_authority: http://dx.doi.org
identifier_product_doi: 10.5067/TERRA/MODIS_OC.2014.0
keywords: Oceans > Ocean Temperature > Sea Surface Temperature
keywords_vocabulary: NASA Global Change Master Directory (GCMD) Science Keywords
data_bins: 15961783
data_minimum: -1.43500006198883
data_maximum: 35.131519317627
sst01.d<-ncvar_get(sst01, 'sst4'); ### check the variable name is correct, can be 'sst4', or other names
sst01.d[sst01.d>42.00072]<-NA ### beware this value should be changed to NA in the matrix
sstlon <- ncvar_get(sst01, 'lon')
sstlat <- ncvar_get(sst01, 'lat')
An image with a colour pallette can be generated using imagep from the oce package
imagep(sstlon, sstlat, sst01.d, col=oceColorsTemperature(255), filledContour=T, missingColor=0, zlab='sst', xlim=c(140,165), ylim=c(-46, -25), zlim=c(10,28), asp=1);
map('worldHires', xlim=c(145,160), ylim=c(-46, -20), fill=T, col="lightgrey", add=T, border=NA)
title(xlab = 'Longitude',
ylab = 'Latitude')
Next zoom in on the RV Investigator sites, and plot them on the map. BE patient.
imagep(sstlon, sstlat, sst01.d, col=oceColorsTemperature(255), filledContour=T, missingColor=0, zlab='sst (˚C)', xlim=c(145,160), ylim=c(-46, -25), zlim=c(10,28), asp=1)
points(in16.sites$long, in16.sites$lat, pch=16, col=in16.sites$col, cex=2)
map('worldHires', xlim=c(145,160), ylim=c(-46, -20), fill=T, col="lightgrey", add=T, border=NA)
title(xlab = 'Longitude',
ylab = 'Latitude')
That is it. Well Done!
Save the file. Preview the .html version. Keep this tutorial somewhere safe for future reference.